#' 中介MR分析（两步法）
#'
#' @param EO_res 暴露到结局的res
#' @param EM_res 暴露到中介的res
#' @param MO_res 中介到结局的res
#' @param method "1"：delta方法；"2"：误差传播法；"3"：另一篇文献的delta方法；
#'
#' @return 中介效应
#' @export
#'
#' @examples
#' # 中介效应是两个单变量计算的
#' # INDIRECT = TOTAL (exposure -> mediator) x TOTAL (mediator -> outcome)
#'
#' \dontrun{
#' library(Oneclick)
#' library(TwoSampleMR)
#' packageVersion("Oneclick")
#' sessionInfo()
#' set.seed(123)
#'
#' ao<-U1_get_ao()
#'
#' # 暴露，中介，结局的数据库两两互不重叠
#' # 此功能可用于多暴露，多中介，多结局
#'
#' # 暴露
#' exposure<- c('finn-b-HYPOTHYROIDISM')
#'
#'
#' # 中介
#'
#' path <- "F:/OneclickDatabase/Oneclick_GWAS_ID/1400种代谢物"
#' #
#' mediation<- file.path( path ,dir(path) )
#'
#' mediation <- "F:/OneclickDatabase/Oneclick_GWAS_ID/1400种代谢物/GCST90200053.parquet"
#'
#' # 结局
#'
#' outcome <- 'ebi-a-GCST90001390'
#'
#' #### 创建文件夹########################################################
#'
#' dir.create("R.file")
#'
#'
#' #### 总效应 暴露到结局########################################################
#'
#' exp_IVs<- U2_extract_instruments(exposure)
#' save( exp_IVs , file = "R.file/exp_IVs.RData" )
#'
#' Outs<- U3_extract_outcomes_data(outcome = outcome ,exposure_iv = exp_IVs )
#'
#' EO_dat<- U4_harmonise_data(exp_IVs,Outs)
#' save( EO_dat , file = "R.file/EO_dat.RData" )
#'
#' EO_res<-U5_mr(EO_dat,run_mr_presso = FALSE,workers = 1)
#' save( EO_res , file = "R.file/EO_res.RData" )
#'
#'
#'
#' #### 第一步 暴露到中介########################################################
#'
#' Outs <- U3_extract_outcomes_data(outcome = mediation ,exposure_iv = exp_IVs )
#'
#' EM_dat<- U4_harmonise_data(exp_IVs,Outs)
#' save( EM_dat , file = "R.file/EM_dat.RData" )
#'
#' EM_res<-U5_mr(EM_dat,run_mr_presso = FALSE,workers = 1)
#' save( EM_res , file = "R.file/EM_res.RData" )
#'
#'
#'
#' #### 第二步 中介到结局########################################################
#'
#' med_IVs<- U2_extract_instruments(mediation,p = 1e-5,clump = "local")
#' save( med_IVs , file = "R.file/med_IVs.RData" )
#'
#'
#' Outs <- U3_extract_outcomes_data(outcome = outcome ,exposure_iv = med_IVs )
#'
#' MO_dat<- U4_harmonise_data(med_IVs,Outs)
#' save( MO_dat , file = "R.file/MO_dat.RData" )
#'
#' MO_res<-U5_mr(MO_dat,run_mr_presso = FALSE,workers = 1)
#' save( MO_res , file = "R.file/MO_res.RData" )
#'
#'
#' #### 中介效应分析 ########################################################
#'
#' mediation_effect <- M8_two_step(EO_res,EM_res,MO_res)
#'
#' save( mediation_effect , file = "R.file/mediation_effect.RData" )
#'
#' colnames(mediation_effect)
#' sig<-subset(mediation_effect, EO_pval<0.05 & EM_pval<0.05 & MO_pval<0.05
#'             & indirect_pval<0.05  )
#'
#'
#'
#'
#'
#' }
#'



M8_two_step<-function(EO_res,EM_res,MO_res,method="1"){

  #暴露到结局
  EO_res$EO_b <- ifelse(!is.na(EO_res$`b_Inverse variance weighted`),
                        EO_res$`b_Inverse variance weighted`,
                        EO_res$`b_Wald ratio`)
  EO_res$EO_se <- ifelse(!is.na(EO_res$`se_Inverse variance weighted`),
                         EO_res$`se_Inverse variance weighted`,
                         EO_res$`se_Wald ratio`)
  EO_res$EO_pval <- ifelse(!is.na(EO_res$`pval_Inverse variance weighted`),
                           EO_res$`pval_Inverse variance weighted`,
                           EO_res$`pval_Wald ratio`)


  #暴露到中介
  EM_res$EM_b <- ifelse(!is.na(EM_res$`b_Inverse variance weighted`),
                        EM_res$`b_Inverse variance weighted`,
                        EM_res$`b_Wald ratio`)
  EM_res$EM_se <- ifelse(!is.na(EM_res$`se_Inverse variance weighted`),
                         EM_res$`se_Inverse variance weighted`,
                         EM_res$`se_Wald ratio`)
  EM_res$EM_pval <- ifelse(!is.na(EM_res$`pval_Inverse variance weighted`),
                           EM_res$`pval_Inverse variance weighted`,
                           EM_res$`pval_Wald ratio`)
  EM_res<- dplyr::rename( EM_res,id.mediation = id.outcome )
  EM_res<- dplyr::rename( EM_res,mediation = outcome )

  #中介到结局
  MO_res$MO_b <- ifelse(!is.na(MO_res$`b_Inverse variance weighted`),
                        MO_res$`b_Inverse variance weighted`,
                        MO_res$`b_Wald ratio`)
  MO_res$MO_se <- ifelse(!is.na(MO_res$`se_Inverse variance weighted`),
                         MO_res$`se_Inverse variance weighted`,
                         MO_res$`se_Wald ratio`)
  MO_res$MO_pval <- ifelse(!is.na(MO_res$`pval_Inverse variance weighted`),
                           MO_res$`pval_Inverse variance weighted`,
                           MO_res$`pval_Wald ratio`)

  MO_res<- dplyr::rename( MO_res,id.mediation = id.exposure )
  MO_res<- dplyr::rename( MO_res,mediation = exposure )

  merge_df <- merge(EO_res, EM_res, by=c("id.exposure","exposure"))
  merge_df <- dplyr::select(merge_df,c("id.exposure","exposure",
                                       "id.mediation","mediation",
                                       "id.outcome","outcome",
                                       'EO_b','EO_se','EO_pval',
                                       'EM_b','EM_se','EM_pval') )

  merge_df <- merge(merge_df, MO_res, by=c("id.mediation","mediation",
                                           "id.outcome","outcome"))

  merge_df <- dplyr::select(merge_df,c("id.exposure","exposure",
                                       "id.mediation","mediation",
                                       "id.outcome","outcome",
                                       'EO_b','EO_se','EO_pval',
                                       'EM_b','EM_se','EM_pval',
                                       'MO_b','MO_se','MO_pval') )


  pb <- progress::progress_bar$new(total = nrow(merge_df))
  mediation_effect_all<-c()
  for (i in 1:nrow(merge_df) ) {

    if(method=="1" ){
    mediation_effect <- tryCatch(
      with(merge_df, mediation_prop(EM_beta=EM_b[i],EM_se=EM_se[i],
                                    MO_beta=MO_b[i], MO_se= MO_se[i],
                                    EO_beta=EO_b[i], EO_se= EO_se[i] ) ),
      error = function(e) return(NA)
    )}

    if(method=="2" ){
      mediation_effect <- tryCatch(
        with(merge_df, product_method_PoE(EM_beta=EM_b[i],EM_se=EM_se[i],
                                      MO_beta=MO_b[i], MO_se= MO_se[i] ) ),
        error = function(e) return(NA)
      )}

    if(method=="3" ){
      mediation_effect <- tryCatch(
        with(merge_df, product_method_Delta(EM_beta=EM_b[i],EM_se=EM_se[i],
                                          MO_beta=MO_b[i], MO_se= MO_se[i] ) ),
        error = function(e) return(NA)
      )}


    mediation_effect$id.exposure <- merge_df$id.exposure[i]
    mediation_effect$id.mediation <- merge_df$id.mediation[i]
    mediation_effect$id.outcome <- merge_df$id.outcome[i]
    mediation_effect$exposure <- merge_df$exposure[i]
    mediation_effect$mediation <- merge_df$mediation[i]
    mediation_effect$outcome <- merge_df$outcome[i]

    mediation_effect_all<-plyr::rbind.fill(mediation_effect_all,mediation_effect )

    pb$tick()
  }

  merge_df <- merge(merge_df, mediation_effect_all, by=c("id.exposure","exposure",
                                                         "id.mediation","mediation",
                                                         "id.outcome","outcome"))
  return(merge_df)


}


mediation_prop <- function(EM_beta, EM_se, MO_beta, MO_se, EO_beta, EO_se ) {

  indirect_beta <- EM_beta*MO_beta
  m1 <- eval(D(expression(EM_beta*MO_beta), "EM_beta"))
  m2 <- eval(D(expression(EM_beta*MO_beta), "MO_beta"))
  indirect_se <- sqrt((m1^2)*EM_se^2 + (m2^2)*MO_se^2)

  prop_mediated_beta <- indirect_beta / EO_beta
  m3 <- eval(D(expression(indirect_beta / EO_beta), "indirect_beta"))
  m4 <- eval(D(expression(indirect_beta / EO_beta), "EO_beta"))
  prop_mediated_se <- sqrt((m3^2)*indirect_se^2 + (m4^2)*EO_se^2)

  res <- data.frame(
    indirect_beta = indirect_beta,
    indirect_se = indirect_se)

  res$indirect_loci = res$indirect_beta-1.96*res$indirect_se
  res$indirect_upci =  res$indirect_beta+1.96* res$indirect_se
  res$indirect_or = exp( res$indirect_beta)
  res$indirect_or_loci = exp( res$indirect_loci)
  res$indirect_or_upci = exp( res$indirect_upci)
  res$indirect_pval = as.numeric(2 * stats::pnorm(abs( res$indirect_beta/ res$indirect_se),lower.tail=FALSE))
  res$prop_mediated_beta =  prop_mediated_beta
  res$prop_mediated_se = prop_mediated_se
  res$prop_mediated_loci = res$prop_mediated_beta - 1.96*res$prop_mediated_se
  res$prop_mediated_upci = res$prop_mediated_beta + 1.96*res$prop_mediated_se
  res$prop_mediated_pval = as.numeric(2 * stats::pnorm(abs(res$prop_mediated_beta/res$prop_mediated_se),lower.tail=FALSE))

  message("使用delta方法计算中介效应.(Doob JL. The limiting distributions of certain statistics. Ann Math Stat. 1935;6(3):160–9),
代码原文PMID: 35104295; delta方法可以在代码原文的第62篇引文查看")


  return (res)
}


product_method_PoE <- function(EM_beta, EM_se, MO_beta, MO_se, verbose=F){
  # calculate indirect effect of exposure on outcome (via mediator)
  # i.e. how much mediator accounts for total effect of exposure on outcome effect

  # this function can run either of two method, depending on what MO data df you supply

  # method 1
  # INDIRECT = TOTAL (exposure -> mediator) x TOTAL (mediator -> outcome)
  # method 2
  # INDIRECT = TOTAL (exposure -> mediator) x DIRECT (of mediator , mvmr)


  # calculate indirect effect beta
  EO_beta <- EM_beta * MO_beta

  if (verbose) {print(paste("Indirect effect = ", round(EM_beta, 2)," x ", round(MO_beta,2), " = ", round(EO, 3)))}


  # calculate SE of indirect effect
  ### using propagation of errors method
  # SE of INDIRECT effect (difference) = sqrt(SE EM^2 + SE MO^2)
  EO_se = round(sqrt(EM_se^2 + MO_se^2), 4)
  if (verbose) {print(paste("SE of indirect effect = sqrt(",
                            round(EM_se, 2),"^2 + ", round(MO_se,2),
                            "^2) = ", indirect_se))}


  # put data into a tidy df
  df <-data.frame(b= EO_beta,
                  se = EO_se)

  df$pval<- as.numeric(2 * stats::pnorm(abs(df$b/df$se),lower.tail=FALSE))

  # calculate CIs and OR
  df$lo_ci    <- df$b - 1.96 * df$se
  df$up_ci    <- df$b + 1.96 * df$se
  df$or        <-  exp(df$b)
  df$or_lci95 <- exp(df$lo_ci)
  df$or_uci95 <- exp(df$up_ci)



  message("使用误差传播法计算中介效应,使用误差传播法计算中介效应的代码原文https://www.nature.com/articles/s42003-022-03272-5")

  return(df)
}


product_method_Delta <- function(EM_beta, EM_se, MO_beta, MO_se, verbose=F){
  # calculate indirect effect of exposure on outcome (via mediator)
  # i.e. how much mediator accounts for total effect of exposure on outcome effect

  # this function can run either of two method, depending on what MO data df you supply

  # method 1
  # INDIRECT = TOTAL (exposure -> mediator) x TOTAL (mediator -> outcome)
  # method 2
  # INDIRECT = TOTAL (exposure -> mediator) x DIRECT (of mediator , mvmr)

  ins <- rownames( installed.packages() )

  if(!"RMediation" %in% ins ){
    install.packages("RMediation")
  }

  # calculate indirect effect beta
  EO <- EM_beta * MO_beta

  if (verbose) {print(paste("Indirect effect = ", round(EM_beta, 2)," x ", round(MO_beta,2), " = ", round(EO, 3)))}


  # Calculate CIs using RMediation package
  CIs = RMediation::medci(EM_beta, MO_beta, EM_se, MO_se, type="dop")

  # put data into a tidy df
  df <-data.frame(b = EO,
                  se = CIs$SE,
                  lo_ci = CIs[["95% CI"]][1],
                  up_ci= CIs[["95% CI"]][2])

  df$pval<- as.numeric(2 * stats::pnorm(abs(df$b/df$se),lower.tail=FALSE))
  # calculate OR
  df$or        <-  exp(df$b)
  df$or_lci95 <- exp(df$lo_ci)
  df$or_uci95 <- exp(df$up_ci)

  message("使用delta法计算中介效应,使用delta法计算中介效应的代码原文https://www.nature.com/articles/s42003-022-03272-5,
delta法的R包的原文https://link.springer.com/article/10.3758/s13428-011-0076-x")

  return(df)
}


